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POINT PROCESS MODELING OF WILDFIRE HAZARD 
IN LOS ANGELES COUNTY, CALIFORNIA 
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The Burning Index (BI) produced daily by the United States 
government's National Fire Danger Rating System is commonly used 
in forecasting the hazard of wildfire activity in the United States. 
However, recent evaluations have shown the BI to be less effective 
at predicting wildfires in Los Angeles County, compared to simple 
point process models incorporating similar meteorological informa- 
tion. Here, we explore the forecasting power of a suite of more com- 
plex point process models that use seasonal wildfire trends, daily 
and lagged weather variables, and historical spatial burn patterns as 
covariates, and that interpolate the records from different weather 
stations. Results are compared with models using only the BI. The 
performance of each model is compared by Akaike Information Cri- 
terion (AIC), as well as by the power in predicting wildfires in the 
historical data set and residual analysis. We find that multiplicative 
models that directly use weather variables offer substantial improve- 
ment in fit compared to models using only the BI, and, in particular, 
models where a distinct spatial bandwidth parameter is estimated for 
each weather station appear to offer substantially improved fit. 

1. Introduction. This paper explores the use of space-time point pro- 
cess models for the short-term forecasting of wildfire hazard in Los Angeles 
County, California. The region is especially well suited to such an analysis, 
since the Los Angeles County Fire Department and Department of Public 
Works have collected and compiled detailed records on the locations burned 
by large wildfires dating back over a century. The landscape in Los Angeles 
County is uniquely vulnerable to high intensity crown-fires, largely because 
the predominant local vegetation consists of dense, highly flammable con- 
tiguous chaparral shrub [Keeley (2000)]. In addition, the dry summers and 
early autumns in Los Angeles County are typically followed by high winds 



Received April 2009; revised August 2010. 

Key words and phrases. Burning index, conditional intensity, point process, residual 
analysis. 

This is an electronic reprint of the original article published by the 

Institute of Mathematical Statistics in The Annals of Applied Statistics. 

2011, Vol. 5, No. 2A, 684-704. This reprint differs from the original in pagination 

and typographic detail. 



1 



2 



H. XU AND F. P. SCHOENBERG 



known locally as Santa Ana winds [Keeley and Fotheringham (2003)]. These 
offshore winds reach speeds exceeding 100 kph at a relative humidity below 
10%, and are annual events lasting several days to several weeks, creating 
the most severe fire weather in the United States [Schroeder et al. (1964)]. 

In order to forecast wildfire hazard, the United States National Fire Dan- 
ger Rating System (NFDRS), created in 1972, produces several daily indices 
that are designed to aid in planning fire control activities on a fire protec- 
tion unit [Deeming et al. (1977); Bradshaw et al. (1983); Burgan (1988)]. 
These include the Occurrence Index, the Burning Index (BI), and the Fire 
Load Index. These indices are derived from three fire behavior components — 
a Spread Component, an Energy Release Component, and an Ignition Com- 
ponent, that are in turn computed based on fuel age, environmental pa- 
rameters (slope, vegetation type, etc.), and meteorological variables such as 
wind, temperature, and relative humidity. Local wildfire management agen- 
cies may combine these components in different ways or calibrate the in- 
herent parameters to adapt the system to the local environment for wildfire 
hazard assessment. 

Fire managers use this information in making decisions about the appro- 
priateness of prescribed burning or alerts for increased preparedness, both 
in terms of fire suppression staffing and fire prevention activities. Since fire- 
line intensity is an important factor in predicting fire containment and the 
likelihood of fire escape, the Burning Index is the rating of most interest to 
many fire managers [Schoenberg et al. (2010)]. This is especially the case for 
natural crown-fire ecosystems such as southern California shrublands, where 
BI is commonly employed to assess fire danger [Mees and Chase (1991)]. In- 
deed, in Los Angeles County, as well as at least 90% of counties nationwide, 
the BI is the index primarily used by fire department officials as a measure 
of overall wildfire hazard, and its use has been justified largely based on its 
observed empirical correlation with wildfire incidence and burn area in dif- 
ferent regions [Haines et al. (1983); Haines, Main and Simard (1986); Mees 
and Chase (1991); Andrews, Loftsgaarden and Bradshaw (2003)]. However, 
several recent investigations have shown that the BI is far from an ideal 
predictor of wildfire incidence in Los Angeles County; Schoenberg et al. 
(2010) showed that a simple point process model, which used only the same 
weather variables as those incorporated by the BI, vastly outperformed the 
BI in terms of predictive efficacy in Los Angeles County, using historical data 
from 1977-2000. In fact, the simple model in Schoenberg et al. (2010) not 
only offered improvement in terms of likelihood scores such as the Akaike 
Information Criterion (AIC), but the study suggested that substantial im- 
provement in short-term forecasting could be achieved by the simple model 
using the weather variables directly, compared to a point process model that 
interpolates BI measurements. 



POINT PROCESS MODELS OF WILDFIRES 



3 



Here, we adopt the same basic modeling framework of Schoenberg et al. 
(2010), but extend the models in two important ways. First, we consider not 
only daily weather variables but also additional covariates with management 
relevance, such as historical spatial burn patterns and wind direction, using 
the directional kernel regression method described in Schoenberg and Xu 
(2008). Second, unlike the simple models of Mees and Chase (1991) and 
Schoenberg et al. (2010) that average daily weather variables over weather 
stations within Los Angeles County, here we explore models that interpolate 
the records from different weather stations, weighting these data based on 
their spatial distance from the location where wildfire hazard is to be esti- 
mated. Thus, the models considered here should have more direct relevance 
for forecasting wildfire hazard in precise spatial locations within Los An- 
geles County, compared to previous work that essentially averaged weather 
variables and hazard estimates over Los Angeles County as a whole. As with 
Schoenberg et al. (2010), our results are compared with models using the 
BI measurements recorded at each of the weather stations, so that the ef- 
fectiveness of the BI in summarizing the wildfire hazard as a function of the 
weather variables may be assessed. 

While alternative models may be more useful for forecasting long-term 
wildfire hazard, that is, estimating the number of wildfires occurring within 
a month, season, or year, the focus here is on forecasting short-term wild- 
fire hazard, that is, the probability of a wildfire occurring within a specific 
day. To compare the overall performance of the models considered, we em- 
ploy diagnostics including likelihood-based numerical summaries such as the 
Akaike Information Criterion (AIC), as well as power diagrams summariz- 
ing the predictive efficacy of each model for short-term forecasting. Residual 
analysis is also used to highlight specific areas and times where the perfor- 
mance of a model is poor and to suggest areas for improvement. 

The paper proceeds as follows. Section 2 describes the wildfire and weather 
data that are used in the analysis. The models used, as well as methods for 
their estimation, are outlined in Section 3, and methods for goodness-of-fit 
assessment are discussed in Section 4. Section 5 presents the main results, 
and a discussion is given in Section 6. 

2. Data. 

2.1. Wildfire data. Los Angeles County is an ideal test site for mod- 
els for wildfire hazard, with detailed wildfire data having been collected 
and compiled by various agencies, including the Los Angeles County Fire 
Department (LACFD) and the Los Angeles County Department of Public 
Works, the Santa Monica Mountains Recreation Area, and the California 
Department of Forestry and Fire Protection. Regional records of the occur- 
rence of wildfires date back to 1878, and include information on each fire, 
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including its origin date, the polygonal outline of the resulting area burned, 
and the centroidal location of this polygon. LACFD officials have noted that 
the records prior to f950 are believed to be complete for fires greater than 
0.405 km 2 (100 acres), and data since 1950 are believed to be complete 
for fires burning greater than 0.0405 km 2 , or 10 acres [Schoenberg et al. 
(2003)]. As in Schoenberg et al. (2010), our analysis in this paper is focused 
primarily on models for the occurrences of the 592 wildfires burning at least 
0.0405 km 2 recorded between January 1976 and December 2000. The daily 
burn area is highly right-skewed and closely follows the tapered Pareto dis- 
tribution [Schoenberg, Peng and Woods (2003)]. For further details, images 
of the spatial locations of these wildfires, and information about missing 
data, see Peng, Schoenberg and Woods (2005). 

2.2. Meteorological data. Since 1976, daily meteorological observations 
from the Remote Automatic Weather Stations (RAWS) were archived across 
the United States. The analysis here is based on sixteen RAWS located 
within Los Angeles County, California. The RAWS record daily measures of 
many meteorological variables, including air temperature, relative humidity, 
precipitation, wind speed, and wind direction [Warren and Vance (1981)]. 
Summaries of these records are collected daily at 1300 hr and transmitted 
by satellite to a central archiving station. These daily RAWS data are used 
as inputs by the NFDRS in order to construct fire behavior components 
that are in turn combined to construct the BI. It should be noted that data 
were missing on certain days for several of the 16 RAWS, though the biases 
resulting from such missing data are likely to be small; see Peng, Schoenberg 
and Woods (2005) for details. 

3. Methodology. We follow previous research including Schoenberg et al. 
(2010) in modeling the catalog of wildfire centroids in Los Angeles County 
as a realization of a point process that may depend on daily meteorological 
variables. We begin with a basic reference model using merely a spatial back- 
ground rate and seasonal component, and a model using the Burning Index 
in addition to the spatial and seasonal background rates. We then introduce 
competing models that use daily meteorological variables recorded at the 
RAWS, and extend the research of Schoenberg et al. (2010) by including 
additional covariates, such as wind direction and fuel age. Further, instead 
of averaging daily weather variables or the Burning Index over all weather 
stations within Los Angeles County, here we explore methods of obtaining 
an estimated spatial intensity at any location x on any particular day by 
interpolating the meteorological variables from different weather stations, 
weighting each record based on its distance from the location x in question. 
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3.1. A review of point process modeling. A spatial-temporal point pro- 
cess -/V is mathematically defined as a random measure on a spatial-temporal 
region S, taking values in the nonnegative integers Z + or infinity [Daley 
and Vere- Jones (2003)]. In this framework the measure N(A) represents the 
number of points falling in the subset A of S. Since any analytical spatial- 
temporal point process is characterized uniquely by its associated condi- 
tional rate (or intensity) A(s), assuming it exists, modeling of such point 
processes is typically performed by specifying a parametric model for this 
rate. For the case where the spatial region is planar, for any point t in time 
and location (x,y) in the plane, the conditional rate is defined as a limiting 
frequency at which events are expected to occur within time range (t, t + At) 
and rectangle (x, x + Ax) x (y,y + Ay), conditional on the prior history, H t , 
of the point process up to time t. For references on space-time point pro- 
cesses and conditional rates, see, for example, Daley and Vere-Jones (1988) 
or Schoenberg, Brillinger and Guttorp (2002). 

Given a parametric function for X(t,x,y), estimates of the parameters 9 
may be obtained by maximizing the log-likelihood function [see Schoenberg, 
Brillinger and Guttorp (2002), page 1576, or Daley and Vere-Jones (2003), 
page 232, equation 7.24]: 

L{9)= [ ' [ [ log[X(t,x,y;9)]dN(t,x,y)- f * / / X(t,x,y; ,9) dydx dt 

JT Jx Jy JT Jx Jy 

= }log\(ti,Xi,yi;9)- / / X(t,x,y;9)dydxdt. 

i = l J T JxJy 

In the case of a Poisson process, the intuition behind this formula is that 
niLi ^(ti,Xi, yf, 9) reflects the likelihood associated with the observed events, 
and exp{ J^T 1 J x J y X(t, x, y; 9) dy dx dt} represents the probability of no events 
in any other portions of the spatial-temporal region, the full likelihood is the 
product of these two terms, and the logarithm of this product yields L(9) 
above. Under rather general conditions, the maximum likelihood estimates 
(MLEs) are consistent, asymptotically normal, and efficient [Ogata (1978)], 
and estimates of their variance can be derived from the negative of the 
diagonal elements of the inverse Hessian of the likelihood function [Ogata 
(1978), Rathbun and Cressie (1994)]. In most cases, explicit solutions for 
MLEs are not available and iterative numerical optimization methods are 
used instead. 



3.2. A simple reference model. In this analysis we explore several spatial- 
temporal point process models for predicting wildfire occurrence rates. As an 
initial baseline model, one may consider an inhomogeneous Poisson process, 
where the conditional intensity at time t and at location (x,y) depends only 
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on the season associated with time t, as well as the background rate m(x,y) 
of wildfires for the location in question. That is, one may consider a baseline 
model such as 



where 7 and a are parameters to be estimated in modeling fitting. 

Parametric or nonparametric methods can be used to estimate the sea- 
sonal pattern S(t) and spatial background m(x,y). While nonparametric 
methods can be especially flexible for estimating complex patterns such as 
spatial burn averages, a possible drawback to such methods is their poten- 
tial for overfitting, particularly when the same data are used for fitting and 
evaluation of the fit of the model. As in Schoenberg et al. (2010), we propose 
estimating the spatial background m(x,y) for fires between 1976-2000 by 
kernel smoothing the centroidal locations of wildfires recorded during the 
previous 25 years, that is, from January 1950 to December 1975. That is, 



where K is a kernel function, j3 m is a bandwidth to be estimated in modeling 
fitting, (xj,yj) indicates the spatial coordinates of the jth wildfire between 
1950 and 1975, no is the number of observed 1951-1975 wildfire occurrences, 
and — (xj,yj)\\ is the Euclidian distance between (x,y) and (xj,yj). 

Standard kernel functions can be used, and attention is usually limited to 
functions that are unimodal, symmetric about zero, and that integrate to 1, 
such as the Gaussian density of the Epanechnikov kernel [Hardle (1994)]. 
It is well known that the results are far more sensitive to the choice of 
bandwidth than the choice of kernel function, and much research has focused 
on automated methods for choosing bandwidth parameters, including cross- 
validation, penalty functions, and plug-in methods [Silverman (1986); Hardle 
(1994)]. Here, since the data (xj,yj) used in the estimation of m(x,y) is 
distinct from that used in the rest of the model fitting and in the evaluation, 
the problem of overfitting is far less severe, and the bandwidth parameter 
may simply be fitted by maximum likelihood. 

Figure 1 shows an estimate of the spatial background rate m(x,y), with 
bandwidth estimated by maximum likelihood. One sees the general pattern 
of fire activity in Los Angeles County during 1951-1975, with most fires 
occurring in the Angeles National Forest, as well as parts of the Los Padres 
National Forest and the Santa Monica Mountains, while many other wildfires 
were located in or near Buckweed, Santa Clarita, and Glendale, California. 

Helmers, Magku and Zitikis (2003) propose a kernel-based estimate for 
the consistent estimation of a seasonal time series. Here, in order to safe- 
guard against overfitting, we propose estimating the seasonal pattern S(t) 




Xi(t,x,y) =jm(x,y) + aS{t) 
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Fig. 1. Spatial background rate m(x,y), with centroid locations of wildfires occurring 
during 1878-1976. (The spatial bandwidth /3 m is 0.6 miles). 



describing the overall seasonal variation of wildfire activity in a fashion simi- 
lar to that used for the spatial background rate, that is, by kernel smoothing 
the times of wildfires during previous years: 

In the above equation, T*(t) represents the date within the year associated 
with time t, that is, T*(t) is the number of days since the beginning of 
the year for time t, tj is the time of the jth wildfire occurrence in the 
data set (1950-1975), and fit is a bandwidth parameter to be estimated. 
A wrapped kernel function K should be used so that, for instance, January 1 
and December 31 are treated as one day apart. The bandwidth may be 
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Fig. 2. Estimate of the seasonal pattern S(t) for model (1), using kernel regression with 
a wrapped Gaussian kernel and bandwidth fit = 9.86 days, estimated by MLE. 



estimated by maximum likelihood, fitting the kernel smoothing of the 1950- 
1975 data to the 1976-2000 data set. This procedure may be preferable for 
relatively small data sets such as the one considered here in order to prevent 
overfitting. 

Figure 2 displays the smoothed function Sit) applied to wildfire incidence 
in Los Angeles County from January 1950 to December 1975, with band- 
width estimated by MLE by fitting the resulting function to wildfire data 
from 1976-2000. It is evident that the mean number of wildfires is highest 
between July and October and rapidly decreases during November and De- 
cember, reaching its minimum in January and February. Schoenberg et al. 
(2010) pointed out that the Burning Index typically assumes moderate val- 
ues in December, January, and February, though few wildfires occur during 
these months. 

3.3. A point process model using Burning Index. To evaluate the po- 
tential of the Burning Index (BI) in predicting wildfire incidence, one may 
consider a model such as 

(2) \ 2 (t,x,y)=jm(x,y) + aS(t) + /j, B iB(t,x,y) 

for some function B(t,x,y) which interpolates the BI records at time t and 
location (x,y), since BI records are only available at fixed RAWS sites. 
Different methods of interpolation are possible. One possibility is to average 
the BI records on day t, weighing each by the distance between the RAWS 
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and the location (x,y) in question. That is, 




\\(x,y) - (a; s ,y s )|| 



) 



BI(t,s) 



where BI(i, s) is the BI value recorded at time t from the sth station, (x s ,y s ) 
are the coordinates of the sth station, St represents the collection of stations 
for which BI records are available on day t, and Cbi is a normalizing constant 
given by 



3.4. Models using spatial interpolation of meteorological variables, includ- 
ing wind speed and wind direction. As an alternative to the model (2) incor- 
porating BI measurements, one may instead consider examining the direct 
impact on wildfire hazard estimates of meteorological variables used in the 
computation of the BI, by replacing the function B(t,x,y) in (2) by func- 
tions of the meteorological variables themselves. That is, one may consider 
models such as 



where Fi(t,x,y) takes into account the contribution of temperature (T), 
relative humidity (H), wind speed (W), and precipitation (P) at time t from 
each RAWS where the data are available. 

Since nonlinearities have been detected in the dependence of burn area on 
climatic variables [Schoenberg et al. (2003)], one may wish to avoid simple 
averaging of the meteorological variables in estimating wildfire hazard. In- 
stead, one option is to describe the association between each climatic variable 
and wildfire burn area by an explicit function g and weight the information 
from each RAWS by the distance to the point (x,y) to be estimated using 
kernel smoothing. This suggests a model such as 




(3) 



h{t,x,y) ='ym(x,y) + aS{t) + F 1 {t,x,y) 
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where T(t,s), H(t,s), W(t,s), and P(t,s) are records of temperature, rela- 
tive humidity, directed wind speed, and precipitation, respectively, on day t 
at the sth RAWS, the parameters fix, fJ>R, f^w, and \ip represent weights as- 
sociated with these meteorological variables, fix, Ph, Pw, and f3p are band- 
widths to be estimated, and Cx, Ch,Cw, and Cp are normalizing constants. 

Note that the bandwidth parameters are somewhat different here than in 
ordinary kernel regression models. While ordinarily in kernel regression or 
kernel density estimation bandwidth parameters may not typically be esti- 
mated by maximum likelihood because the likelihood would tend to increase 
as the bandwidth shrinks to [Silverman (1986)], here this is not the case. 
Instead, the bandwidth parameters fix, flu, Pw , and l3p in model (3) merely 
control the spheres of influence of the relative weather stations in terms of 
the impact of each on wildfire hazard. That is, if fix is small, for instance, 
then each RAWS station's recorded temperature will affect the wildfire in- 
cidence more locally, whereas if (3x is very large, then the wildfire hazard at 
any particular location will depend more closely on the average temperature 
throughout Los Angeles County. 

Functional forms can be suggested for gT,9H,9w, and gp, by individually 
examining the empirical relationship between daily area burned and each 
of these variables. In order to smooth these relationships, one possibility 
would be to use local linear regression or segmented regression, since the 
relationships between wildfire burn area and temperature, precipitation, and 
other weather variables appear to have thresholds [Schoenberg et al. (2003)]. 
Another possibility is to use kernel regression of daily area burned on the 
average temperature, relative humidity, and precipitation over all RAWS, 
respectively. For instance, the impact of temperature may be estimated via 



where Aj is the area burned on the jth day during 1976-2000, Tj is the aver- 
age temperature readings over all RAWS on that day, hx is the bandwidth 
of the kernel regression which can be selected by methods such as cross- 
validation or the plug-in method [Silverman (1986)], and n\ is the number 
of days with records during this period. Figure 3 displays such kernel regres- 
sion estimates of gx, gu, and gp. Not surprisingly, one sees that daily area 
burned generally increases as temperature increases, and decreases as rela- 
tive humidity and precipitation increase, though some local fluctuations are 
seen in the kernel regressions on temperature and relative humidity. These 
fluctuations are likely attributable to the high variability of the estimates 
due to the relatively small sample of large fires contained in the catalog. 

Special care should be taken in estimating g\y, since wind is directional, 
and this direction may provide important information related to wildfire 
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Fig. 3. Kernel regression estimates of the relationship between daily burn area and 
(a) temperature, (b) relative humidity, and (c) precipitation. Gaussian kernels are used, 
bandwidths are estimated by cross-validation, and edge correction is performed via reflec- 
tion [Silverman (1986)]. 



incidence. One possible way to estimate the relationship between daily area 
burned and directional wind speed is via directional kernel regression, as 
outlined in Schoenberg and Xu (2008). An example of a two-dimensional 
directional kernel is the von Mises distribution suggested by Mardia and 
Jupp (2000): 

«M(0; /*,«) = — l — e — (^), 
2tt1o{k) 

where Iq denotes the modified Bessel function of the first kind and order 0, 
fi is the directional center, and k is known as the concentration parameter. 
Following Schoenberg and Xu (2008), the corresponding two-dimensional 
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Fig. 4. Two-dimensional kernel regression of burn area versus wind speed and wind 
direction. Wind speed and wind direction are represented in a polar system: wind speed 
is represented by the distance to the center and wind direction is represented by the angle. 
The grey scale represents the smoothed average wildfire burn area during 1976-2000. 



kernel regression function gw would then be estimated via 

(W = ^T=ii K (\ W ~ W,\lh w )vM(B - Oj-M, Ko)Aj} 

9W[ ' ' T^LAKdw-Wji/hw^Mie-ej-^Ko)} ' 

where Wj and 9j represent the mean wind speed and wind direction, re- 
spectively, on day j. Cross-validation can be used to optimize the estimates 
of h sp , no, and kq. 

Figure 4 displays a kernel regression estimate of the relationship between 
daily burning area and daily mean wind direction, weighted by wind speed, 
averaged over all 16 RAWS stations. The sharp increase in mean wildfire 
burn area, indicated by darker shading in Figure 4, is very strongly asso- 
ciated with higher wind speeds. In addition, one sees from Figure 4 the 
extent to which winds from the northeast, which are often warm, dry Santa 
Ana winds, are associated with higher burn areas. Since the impact on av- 
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erage wildfire burn area of wind direction might be different at distinct 
weather stations, one might wish to estimate 16 distinct kernel regression 

(s) 

functions g\y , one for each RAWS station s. 

The model (3) described above is additive in each of the weather variables, 
implying that an extreme value in only one weather variable may lead to 
a high estimate of wildfire hazard on the corresponding day, which might 
be questionable. For instance, one would expect few large wildfires occur 
on days when temperatures are extremely high yet there is some moderate 
amount of precipitation and relative humidity. An alternative approach is 
to use a multiplicative component instead, where once again 



3.5. Models allowing records at different RAWS to have different relation- 
ships with wildfire hazard. In both model (3) and model (4), gT,gu, and gp 
may be estimated using kernel regression of burn area on average temper- 
ature, relative humidity, and precipitation over all weather stations, where 
each station has the same regression function. However, because of differ- 
ences in the locations of the weather stations, including differing altitudes 
of these stations, some stations may have lower average temperatures or 
higher relative humidities than others throughout the year. Hence, a partic- 
ular temperature and relative humidity at one station might indicate a very 
different wildfire hazard than the same values observed at a different RAWS 
station. In order to deal with this, one may estimate distinct kernel regres- 
sion curves for each RAWS in model (3) and model (4). That is, one might 
consider 



(4) \ 4 (t,x,y) =jm(x,y) + aS(t) + fiF 2 (t,x,y), 

where now the fire weather (F) term has the multiplicative form 




^2^ V P2 



(5) 



h{t,x,y) ='ym(x,y) + aS(t) + F 3 (t,x,y) 



where 



Fs{t,x,y) 
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or 



(6) 



Ae(i, x, y) = jm(x, y) + as(t) + nF 4 (t, x, y) 



where 



F 4 (t,x,y) 



C 



1 




\\(x,y) - (x s ,y s )\\\ ( s) 




The kernel regression functions such as (t,s) for each station s may be 
estimated as in models (3) and (4), that is, by kernel regression of the total 
daily burn area in Los Angeles County against the temperature at station s. 

3.6. Incorporating fuel age. One may further improve the models by 
adding fuel age as a covariate. Fuel age, or its proxy, the time since the 
location's last recorded burn, appears to have a nonlinear, threshold-type 
relationship with burn area [Peng and Schoenberg (2008)]. Indeed, burn 
area appears to increase steadily with fuel age up to ages of approximately 
20-30 years [Peng and Schoenberg (2008)]. This suggests incorporating the 
contribution of fuel age into model (5) by a truncated linear function, that is, 

(7) A 7 (i, x, y) = jm(x, y) + as(t) + F 3 (t, x, y) + \i D mm{D(t, x, y),ip}, 

where D(t,x,y) is the fuel age at the space-time pixel (t,x,y), and where tp 
is an upper truncation time. Fuel age may be incorporated similarly into 
model (6) as well: 

(8) X 8 (t,x,y) = jm(x,y) + aS(t) +{iF 4 (t,x,y) + fi D min{L»(t, x, y), tp}. 

4. Model assessment. Equations (l)-(8) describe eight point process 
models that may be used to predict wildfire hazard at any time and location 
within Los Angeles County. In order to compare the performance of these 
models, one commonly used method is the Akaike Information Criterion 
(AIC), which is defined as — 2L(/3) + 2p, where L(f3) is the logdikelihood 
and p is the number of fitted parameters in the model. Smaller values of 
AIC indicate better fit. The AIC makes a good trade-off between model 
complexity and overfitting by rewarding a higher likelihood while penalizing 
the addition of more parameters [Akaike (1977)]. 

The predictive capacity of competing point process models may also be 
compared by examining the models' performance on the 1976-2000 wildfire 
data, as suggested in Schoenberg et al. (2010). Consider a grid of space-time 
cells, with each cell's center separated by some distance Ad in space and 
a temporal distance At, and let these cells represent locations and times 
where alarms may potentially be issued. For any such space-time point 
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(t,x,y), one may compute the estimated conditional intensity X(t,x,y) for 
a particular model. Consider issuing an alarm if the value of \(t,x,y) is 
above some certain threshold. We say the alarm is successful if a wildfire 
occurs within the cell; otherwise, it is a false alarm. The false positive rate 
of the alarms, defined as the proportion of cells without wildfires where A 
exceeded the alarm threshold, can be compared to the true positive rate, 
that is, the proportion of wildfires occurring in cells where A exceeded the 
alarm threshold, using traditional Receiver Operating Characteristic (ROC) 
curves. Each possible alarm threshold represents a single point on the ROC 
curve, and the resulting curve summarizes the potential efficacy of a model 
in forecasting wildfires. 

While numerical likelihood scores such as AIC and ROC curves can be 
useful in evaluating the overall performance of a point process model, neither 
method is useful at identifying particular times and locations where a model 
fits poorly or suggesting ways in which a model might be improved. For these 
purposes, it is useful to inspect plots of residuals, which may be defined as the 
difference between the number of events occurring in a certain space-time 
interval and the integral of the estimated conditional intensity over the same 
interval [Baddeley et al. (2005)]. Negative residuals indicate overestimates 
of wildfire hazard, and very large residuals indicate places and times where 
the model underestimated wildfire hazard. 

5. Results. The maximum likelihood estimates of the parameters for the 
models (1)— (8) are listed in Tables 1 and 2. In Tables 1 and 2, the parame- 
ter ip was fixed at 22 years for models 7 and 8, based on Peng and Schoenberg 
(2008); this parameter was also fit by maximum likelihood, yielding very sim- 
ilar results, so, for simplicity, here we report the fit of the model with ip fixed 
at 22 years. The bandwidths in spatial background f3 m range from 0.25 km 
to 1.20 km and the bandwidth in the seasonal component fall within 8.6 to 
34.1 days. The bandwidths related to spatially kernel smoothing the weather 
variables range from 0.024 km to 0.40 km in models (3), (5), and (7), with 
the smallest value for wind speed in model (7) and the largest value corre- 
sponding to relative humidity in model (3). As mentioned in Section 3, these 
bandwidths can perhaps be interpreted as reflecting the scales of influence 
of the weather variables in terms of their effect on wildfire incidence. 

Table 3 presents the relative AIC values for models (l)-(8). For simplicity 
and ease of presentation, the AIC for the best fitting model (8) has been 
subtracted from the AIC of each model. It is evident that the BI model (2) 
offers very substantial improvement over the baseline model (1). However, all 
the other models that use weather information directly have much better fits 
than the BI model (2). The multiplicative model (8) with fuel age appears to 
offer by far the best fit among these models, without using the BI directly, 
and only involves one more fitted parameter than the BI model. 
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Table 1 

Maximum likelihood estimates of scaling parameters 



Model 


7 


a 


Mb 




fJ-H 










(1) 


24 
(1.3) 


2.0 

(0.20) 
















(2) 


4.9 


0.65 


8.6 x 10" 4 
















(0.26) 


(0.039) (2.0 x 10" 5 ) 














(3) 


6.4 


0.66 




0.18 


0.21 


1.0 


0.15 








(0.64) 


(0.043) 




(0.017) 


(0.012) 


(0.071) 


(0.010) 






(1) 


6.4 
(0.31) 


2.7 
(0.24) 












2.1 x 10 3 
(130) 




(5) 


13 




0.60 


0.58 


0.19 


17 


2.1 








(0.84) 




(0.051) 


(0.054) 


(0.011) 


(0.85) 


(0.14) 






(6) 


12 

(0.15) 


1.0 

(0.057) 












5.1 x 10 3 
(260) 




(7) 


6.9 


0.60 




0.19 


0.21 


52 


0.71 




1.0 




(0.44) 


(0.054) 




(0.018) 


(0.020) 


(2.4) 


(0.068) 




(0.066) 


(8) 


4.8 
(0.37) 


0.54 
(0.049) 












1980 x 10 3 
(85 x 10 3 ) 


0.10 
(0.0094) 



All entries have been multiplied by 10 3 for brevity. 



Figure 5 shows a comparison of the predictive efficacy of models described 
in Section 3. Models (6) and (8) vastly outperform the other models. The 
performance was evaluated using a regular space-time grid, so that each 
alarm's success or failure was evaluated over a space-time window with 
Ad = 4.0 km and At = 1.0 day. 



Table 2 

Maximum likelihood estimates of bandwidth parameters 



Model (3 m (km) 


Pt (day) 


/3s (km) 


Ar (km) 


Ph (km) 


[3 W (km) 


(3p (km) 


(3 (km) 


(1) 


1.20 


9.86 
















(0.004) 


(2.8) 














(2) 


0.92 


8.64 


0.40 














(0.00077) 


(3.30) 


(0.0055) 












(3) 


0.36 


29.6 




0.37 


0.40 


0.30 


0.24 






(0.002) 


(5.6) 




(0.002) 


(0.001) 


(0.004) 


(0.001) 




(1) 


0.92 


8.6 












0.03 




(0.1) 


(2.3) 












(0.003) 


(5) 


0.34 


34 




0.31 


0.20 


0.04 


0.28 






(0.002) 


(8.3) 




(0.001) 


(0.001) 


(0.002) 


(0.0005) 




(6) 


0.25 


20 












0.19 




(0.025) 


(9.7) 












(0.003) 


(7) 


0.46 


19 




0.36 


0.39 


0.024 


0.20 






(0.00069) 


(5.2) 




(0.00071) 


(0.0011) 


(3.2e-8) 


(0.00025) 




(») 


0.99 


13 












0.037 




(0.17) 


(2.6) 












(0.00064) 
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Table 3 
Relative AIC values 

Model (1) (2) (3) (4) (5) (6) (7) (8) 

Relative AIC 3783 2859 2509 2631 2223 2209 1308 
p 4 6 12 6 12 6 13 7 



For any given success rate, the models that directly use the meteorological 
data offer substantially fewer false alarms than the model (2) that uses the 
BI. For instance, for a false positive rate fixed at 0.08, model (8) correctly 
signals approximately 29% of the wildfires in the data set, compared to 18% 
for model (2). Model (6), which uses only temperature, relative humidity, 
wind speed, wind direction, and precipitation, but does not use fuel age, 
signals nearly 25% of the wildfires correctly with a false positive rate of 
0.08. Note that this method of evaluating predictive efficacy over a fine 
grid of spatial-temporal locations is rather cumbersome for models (7)-(8), 
due to the need to individually estimate the fuel age associated with each 
wildfire, with respect to each spatial-temporal grid location and time, and 
each such evaluation requires a rather burdensome computation described 
in Peng and Schoenberg (2008). 



2 



1) 

(= 



O 




T" 





model 1 


-2- 


model 2 


-3- 


model 3 




model 4 


-5- 


model 5 


-6- 


model 6 




model 7 


-8- 


model 8 



0.0 



0.2 



0.4 0.6 
False Positive Rate 



0.8 



1.0 



Fig. 5. ROC curves for models (l)-(8), with A t = 1.0 day and Ad = 4.0 km. 
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Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 

Fig. 6. Median absolute value of residuals, by month, using a space-time grid of 
25.6 sqkm x 30.0 days. 

The fit of the models can be evaluated by examining their spatial-temporal 
residuals over a relatively coarse grid. For instance, Figure 6 shows the medi- 
ans of the absolute values of the residuals in each month, for models (l)-(8), 
where each residual is computed over a space-time grid of 25.6 sqkm x 30.0 
days. It is evident that models (3), (5), (6), (7), and (8) outperform the 
other three models, especially in the late Summer and Fall months when 
most wildfires in Los Angeles County occur. The months of October and 
November are especially critical, since Santa Ana winds prevail and can 
cause catastrophic wildfires. Figure 7 shows a spatial plot of the medians 
of the absolute values of the residuals, over all months. From Figure 7 one 
sees that the eight models perform surprisingly comparably in terms of the 
median absolute residual, though models (3)~(8), which use meteorological 
data directly, have better performance than model (2) in the sense that the 
corresponding residuals are generally closer to zero in most areas. These 
residual plots also indicate that several of the models may require improve- 
ment in the Northwest portion of the Angeles National Forest, as well as 
near the border with San Bernardino County on the Eastern side of Los 
Angeles County. From the residuals in Figure 7, one can observe that the 
residuals for models (5), (7), and (8) are highly concentrated around zero 
except a few large values that occur in the Northwest part and east part of 
Los Angeles County. 

6. Discussion. The models explored here use the identical information 
recorded at the RAWS stations and used as inputs into the computation 
of the BI. Hence, it seems relevant to compare the fit of such models with 
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0.0312 
0.026 

0.0208 
I I 0.0156 
I I 0.0104 
I I 0.00521 

n 

Median absolute residual 
fires/km 2 /yr 



Fig. 7. Median absolute value of residuals, by location, using a space-time grid of 
25.6 sqkm x 30.0 days. 

comparable spatial interpolations of BI measurements, and the fact that 
the models using the weather variables directly appear to offer superior fit 
suggests that the BI may not be effective as a short-term forecasting measure 
of wildfire hazard in Los Angeles County. 

It should be noted that the empirical relationship between a fire danger 
rating index such as the BI and wildfire incidence is only one way to eval- 
uate the effectiveness of such an index; alternatives may include assessing 
the cost-effectiveness of staffing or other decisions made based on the index. 
Furthermore, the use of fire danger ratings by fire department officials for 
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wildfire suppression and prevention activities may confound the empirical 
relationship between fire danger ratings and observed wildfire activity. Nev- 
ertheless, most evaluation studies of fire danger rating systems relate such 
indices to ultimate fire responses, including fire incidence and fire size. In- 
deed, Andrews and Bradshaw (1997), whose work was instrumental in the 
current implementation of the BI, suggested that the value of a fire danger 
index be evaluated according to its relationship with fire activity, which may 
be defined as the incidence of large wildfires. Such empirical relationships 
have been used as support for the use of such rating systems for predictive 
purposes [Haines et al. (1983); Haines, Main and Simard (1986); Mees and 
Chase (1991); Mandallaz and Ye (1997a, 1997b); Viegas et al. (1999); An- 
drews, Loftsgaarden and Bradshaw (2003)]. The results here suggest that, 
for the purpose of forecasting wildfire hazard, point process models using 
RAWS records and previous wildfire activity as covariates may represent 
a promising alternative to existing indices that use essentially the same in- 
formation. 

However, we must emphasize that the point process models proposed here 
remain rather simplistic and could potentially be improved by incorporating 
a host of other important variables, such as detailed vegetation type, vege- 
tation cover, soil characteristics, other weather variables such as cloud cover 
and lightning, as well as human factors such as land use and public policy. 
The exclusion of such variables from this analysis is solely motivated by our 
aim to optimize forecasts given the same remote, automatically-recorded in- 
formation used in the computation of the BI. The models considered here 
could also perhaps be improved in various ways. For instance, one might al- 
low long-term temporal trends and /or allow the seasonal component to vary 
from year to year. In addition, one may consider estimating the kernel func- 
tion in models (5) and (6) for each station using only local wildfires close to 
the corresponding station, or perhaps by some more sophisticated weighting 
scheme where nearby fires are given higher weight in the estimation of this 
function. Because daily burn areas are right-skewed [Schoenberg, Peng and 
Woods (2003)], perhaps kernel regressions where the response variable is 
some transformation of the daily burn area might yield superior results. An 
additional important direction for future work is the exploration of similar 
point process models for wildfire occurrences in other locations and for other 
vegetation types or alternative wildfire regimes, as well as the use of such 
models for actual prospective predictions of wildfire activity, rather than 
merely the empirical assessment of goodness of fit to historical data. 
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